%% Construct data for estimating VAR and compute MLE.
YY = data(opt.p+1:end,:); % LHS variables
XX = lagmatrix(data,1:opt.p); % Matrix of regressors
XX = XX(opt.p+1:end,:); % Drop initial missing observations

if opt.const == 1 % Add constant to matrix of regressors if required
    XX = [ones(length(XX),1) XX];   
end

N = size(YY,2); % Number of variables
Ni = length(opt.ivar); % Number of variables of interest
M = N*opt.p + opt.const; % Number of parameters in each equation
T = length(YY); % Number of observations used in estimating VAR

% Compute useful variables and check whether analytical approach of GMM18 
% is applicable.
[opt.calcLRCIR,oneColQ,restr.f,restr.s,convInfo] = initialComp(N,opt,restr);

%% Conduct posterior inference on IRs.
% j-th column of B gives MLE/OLS coefficients of j-th reduced form VAR 
% equation. If constant included, first row corresponds to constant term.
phiHat.B = (XX'*XX)\XX'*YY;
phiHat.S = (YY - XX*phiHat.B)'*(YY-XX*phiHat.B);
phiHat.Sigma = (1/(T-M))*phiHat.S; % Variance matrix of VAR innovations
phiHat.P = (XX'*XX)\eye(M);
phiHat.cholP = chol(phiHat.P,'lower');

% Storage arrays.
rMinPost = zeros(opt.Kdraws,opt.H+1,Ni);
rMaxPost = rMinPost;
rSinglePriorPost = rMinPost;
B = zeros(N*M,opt.Kdraws);
Sigma = zeros(N,N,opt.Kdraws);
convInfo.convInd = zeros(opt.Kdraws,1);

phiDraw = 0; % Counter for no. of draws from posterior of phi
nonStableCount = 0; % Counter for no. of draws with nonstable VAR
nonEmptyDraw = 0; % Counter for no. of draws with non-empty IS

while nonEmptyDraw < opt.Kdraws
    
    % Posterior sampler with noninformative improper (Jeffreys) prior.
    % See the file 'niwPrior.m' for details about implementing a
    % normal-inverse-Wishart prior.
    phi.Sigma = iwishrnd(phiHat.S,T-M);
    phi.Sigmatr = chol(phi.Sigma,'lower');
    phi.Sigmatrinv = phi.Sigmatr\eye(N);
    phi.B = phiHat.B(:) + kron(phi.Sigmatr,phiHat.cholP)*randn(M*N,1);
    phiDraw = phiDraw + 1;

    % Generate coefficients in reduced-form VMA representation and, if 
    % required, long-run (cumulative) IR. Also generate versions of these 
    % objects post-multiplied by Sigmatr, and check invertibility of the 
    % VAR representation.
    [irs,nonStable] = genVMA(phi,opt);
    nonStableCount = nonStableCount + nonStable;
        
    if oneColQ == 1 && opt.whichAlgo == 3
        % Use analytical approach to determine whether IS is nonempty.
        Qempty(phiDraw) = checkEmptyISAnalytical(restr,irs,phi,opt);         
    end
        
    % Use Algorithm 1 (Step 2) to determine whether IS is nonempty and 
    % draw Q0 from space of orthonormal matrices satisfying equality and 
    % sign restrictions. This draw is used to compute the IR under a 
    % single prior for Q.
    [Q0,restr.F,restr.S,QemptyNumerical(phiDraw)] = ...
        drawQ(restr,irs,phi,opt);
    
    %phi.Sigmatr*Q0

    if oneColQ == 0 || opt.whichAlgo ~= 3
        Qempty(phiDraw) = QemptyNumerical(phiDraw);
    end
               
    if Qempty(phiDraw) == 0 % If IS non-empty

    

            % Store reduced-form parameters
            nonEmptyDraw = nonEmptyDraw + 1;
            Sigma(:,:,nonEmptyDraw) = phi.Sigma;
            B(:,nonEmptyDraw) = phi.B;

            % Compute draw from posterior of IRs given single prior for Q.

            for hh = 1:opt.H+1 % For each horizon

                for ii = 1:Ni % For each variable of interest

                    rSinglePriorPost(nonEmptyDraw,hh,ii) = ...
                        irs.vmaGK(opt.ivar(ii),:,hh)*Q0(:,opt.jshock);

                end

            end
            
           % Find upper and lower bounds of IR IS.
            
           if opt.whichAlgo == 1
                
               % Use Algorithm 1 (constrained optimisation) from GK18.
               % Users without the Parallel Computing Toolbox can replace 
               % the function numericalBoundsPar with numericalBounds.
               [rMinPost(nonEmptyDraw,:,:),rMaxPost(nonEmptyDraw,:,:)] = ...
                   numericalBoundsPar(restr,irs,phi,opt,optimOptions);
               
           elseif opt.whichAlgo == 2
               
               % Use Algorithm 2 (a simulation-based method) from GK18.
               [rMinPost(nonEmptyDraw,:,:),rMaxPost(nonEmptyDraw,:,:)] = ...
                   approximateBounds(restr,irs,phi,opt);
                             
           elseif opt.whichAlgo == 3    
               
               % Use analytical results from GMM18.
               [rMinPost(nonEmptyDraw,:,:),rMaxPost(nonEmptyDraw,:,:)] = ...
                   analyticalBounds(restr,irs,phi,opt);

           end
           
           % Check whether IS is convex at each draw of phi using 
           % Proposition 3
           if convInfo.propCase ~= 1 && convInfo.propCase ~= 2
               convInfo.convInd(nonEmptyDraw) = ...
                   checkConvex(restr,phi,opt,convInfo);
           end
           
           if mod(opt.Kdraws-nonEmptyDraw,opt.dispIter) == 0
              
               fprintf('\n%d draws with non-empty identified set remaining...',...
               opt.Kdraws-nonEmptyDraw);
           end     
    end     
end


% Compute posterior mean under single prior.
rSinglePriorMean = permute(mean(rSinglePriorPost,1),[2 3 1]);

% Compute posterior mean bounds.
meanlb = permute(mean(rMinPost,1),[2 3 1]);
meanub = permute(mean(rMaxPost,1),[2 3 1]);
postMeanBoundWidth = meanub - meanlb; % Width of posterior mean bounds.

% Compute robustified credible region for IS.
[credlb,credub] = credibleRegion(rMinPost,rMaxPost,opt);
credWidth = credub - credlb; % Width of robustified credible region

% Compute highest posterior density (HPD) interval under single prior.
[hpdlb,hpdub] = highestPosteriorDensity(rSinglePriorPost,opt);
hpdWidth = hpdub - hpdlb; % Width of highest posterior density regions

% Rows in Table 6 and 7 (left column) for this model:
disp('Table 6')
disp([meanlb(1:6,2) meanub(1:6,2) postMeanBoundWidth(1:6,2) rSinglePriorMean(1:6,2) max(abs(rSinglePriorMean(1:6,2)-meanlb(1:6,2)),abs(meanub(1:6,2)-rSinglePriorMean(1:6,2)))]*100)

disp([meanlb(1:6,3) meanub(1:6,3) postMeanBoundWidth(1:6,3) rSinglePriorMean(1:6,3) max(abs(rSinglePriorMean(1:6,3)-meanlb(1:6,3)),abs(meanub(1:6,3)-rSinglePriorMean(1:6,3)))])

disp([meanlb(1:6,1) meanub(1:6,1) postMeanBoundWidth(1:6,1) rSinglePriorMean(1:6,1) max(abs(rSinglePriorMean(1:6,1)-meanlb(1:6,1)),abs(meanub(1:6,1)-rSinglePriorMean(1:6,1)))]*12)

%disp('Table A.3')
%disp([credlb(1:6,2) credub(1:6,2) credWidth(1:6,2) hpdlb(1:6,2) hpdub(1:6,2) hpdWidth(1:6,2) max(abs(credlb(1:6,2)-hpdlb(1:6,2)),abs(credub(1:6,2)-hpdub(1:6,2)))]*100)
%disp([credlb(1:6,3) credub(1:6,3) credWidth(1:6,3) hpdlb(1:6,3) hpdub(1:6,3) hpdWidth(1:6,3) max(abs(credlb(1:6,3)-hpdlb(1:6,3)),abs(credub(1:6,3)-hpdub(1:6,3)))])
%disp([credlb(1:6,1) credub(1:6,1) credWidth(1:6,1) hpdlb(1:6,1) hpdub(1:6,1) hpdWidth(1:6,1) max(abs(credlb(1:6,1)-hpdlb(1:6,1)),abs(credub(1:6,1)-hpdub(1:6,1)))]*12)



